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Covariant density functional theory, in the framework of self-consistent Relativistic Mean Field 
(RMF) and Relativistic Random Phase approximation (RPA), is for the first time applied to axially 
deformed nuclei. The fully self-consistent RMF-fRRPA equations are posed for the case of axial 
symmetry and non-linear energy functionals, and solved with the help of a new parallel code. Formal 
properties of RPA theory are studied and special care is taken in order to validate the proper 
decoupling of spurious modes and their influence on the physical response. Sample applications to 
the magnetic and electric dipole transitions in ^°Ne are presented and analyzed. 
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I. INTRODUCTION 

New experimental facilities with radioactive nuclear 
beams have stimulated enhanced experimental and the- 
oretical efforts to understand the structure of nuclei, not 
only along the narrow line of stable isotopes, but also in 
areas of large neutron- and proton excess far from the val- 
ley of /3-stability. Beside the investigation of the ground 
state properties of these nuclei, more and more experi- 
mental studies are being devoted to the understanding of 
the properties of excited states in this region. 

On the theoretical side, only very light nuclei can 
be studied in the framework of modern ab initio me- 
thods. Shell model calculations in restricted configu- 
ration spaces provide an accurate description of light 
and medium-heavy nuclei. For the large majority of 
nuclei, however, a quantitative microscopic description 
is only possible using density functional theory (DFT). 
Although DFT can, in principle, provide an exact de- 
scription of the many-body dynamics if the exact density 
functional is known [li , in nuclear physics one is far 
from a microscopic derivation of this functional and in 
addition there is the problem, that in self-bound systems 
density functional theory can only be applied to intrinsic 
densites [1, 01 . The most successful schemes use a phe- 
nomenological ansatz incorporating as many symmetries 
as possible, and adjust the parameters of these function- 
als to ground state properties of characteristic nuclei all 
over the periodic table (for a recent review see Ref. Q). 

Of particular interest are covariant density functionals 
@) 01 because they are based on Lorentz invariance. The 
inclusion of this symmetry not only allows for the descrip- 
tion of the spin-orbit part of the nuclear interaction in a 
natural and consistent way, but it also puts considerable 
restrictions on the number of parameters in the corre- 
sponding functionals, all without reducing the quality of 
the agreement with experimental data. A very success- 
ful example is the Relativistic Hartree-Bogoliubov model 
0, ^], which combines a density dependence through a 
non-linear coupling between the meson fields [al with 
pairing correlations based on an effective interaction of 



finite range Q. 

Excited states are described within this formalism 
by time-dependent density functional theory p^ . In 
the small amplitude limit one obtains the relativistic 
Random Phase Approximation (RRPA) [ll|, El ■ This 
method provides a natural framework to investigate col- 
lective and non-collective excitations of p/i-character. Al- 
though several RRPA irnplementations have been avail- 
able since the eighties [Ij], only very recently RRPA- 
based calculations have reached a level on which a quan- 
titative comparison with experimental data became pos- 
sible JJj] . And even though the self-consistent relativistic 
mean-field (RMF) framework has been employed in many 
studies of deformed nuclei [l^, [11] , applications of the rel- 
ativistic RRPA method have so far restricted to spher- 
ical nuclei. This is also true for non-relativistic density 
functionals, where most of the RPA-calculations are re- 
stricted to spherical nuclei. Only very few deformed RPA 
calculations based on Skyrme (Ojllll or Gogny forces [H 
are available so far. 

On the other hand, it is well known that only semi- 
magic nuclei have spherical shape, and that most of the 
other nuclei in the nuclear chart are deformed. Thus, the 
description of the collective response of these nuclei can 
only be accomplished within a framework where defor- 
mation is explicitly taken into account. From the point 
of view of nuclear structure the motivation for deformed 
RPA calculations is evident. In addition, the nuclear elec- 
tric dipole response obtained in this framework provides 
valuable input for the calculation of important astrophys- 
ical processes [20|], such as the r- or the s-process, that 
pass though large areas of deformed nuclei. 

In this manuscript we report on the extension of rel- 
ativistic RPA theory to axially deformed nuclei and its 
application to the study of collective excitations. In sec- 
tion |lT] we discuss the underlying density functional and 
the derivation of the relativistic RPA equations. Sec- 
tion mil deals with specific aspects of RMF and RPA 
theory in deformed systems and the evaluation of the 
relativistic RPA matrix elements in the basis of axially 
deformed Dirac spinors. Section lTVl is devoted to strength 
functions and sum rules and in section W\ we discuss 
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transition densities in the intrinsic and in the laboratory 
frame. Violations of symmetries and the corresponding 
Goldstone modes are treated in Section IVll and in section 
IVIII we show illustrative applications in ^°Ne, in partic- 
ular its magnetic and electric dipole response. Finally 
section I Villi contains the summary and an outlook. 



II. THE COVARIANT ENERGY DENSITY 
FUNCTIONAL AND THE RELATIVISTIC RPA 
EQUATIONS 

At the moment the most successful density functionals 
in nuclear physics are purely phenomenological. Consid- 
ering from the beginning as many symmetries as possible, 
one starts with a relatively simple ansatz for the energy 
density functional [21, 22, 23], which contains a certain 
number of phenomenological parameters. One then ad- 
justs these parameters to bulk properties of nuclear mat- 
ter and to ground state properties of a few selected finite 
nuclei with spherical shape. These sets are then used over 
the entire nuclear chart. It turns out that for a good de- 
scription of the experiment data, it is crucial to allow for 
a density dependence in this ansatz. the concept of den- 
sity dependence has its origin in more microscopic theo- 
ries of the nuclear many-body system, such as Brueckner 
theory [2^, which leads to a density dependent effective 
interaction in the nuclear interior. In relativistic models 
this density dependence was taken into account in the 
form of a non-linear meson coupling in Ref. Q or in the 
form of density dependent meson-nucleon couplings in 
Ref. [25] . If the ansatz is chosen properly and if the ad- 
justment of the phenomenological parameters is carefully 
done, the quantitative agreement with available experi- 
mental data is remarkable 

In this work we concentrate on relativistic density func- 
tional theory 0]. These functionals are based on 
Lorentz invariance. The basic degrees of freedom are the 
nucleons described by point-like Dirac spinors. In order 
to be consistent with Lorentz invariance and causality, 
one has two possibilities for introducing an interaction 
between these particles. Either one restricts the theory 
to zero range interactions, as it is done in Nambu Jona- 
Lasinio models [1^ , or one allows for the exchange of ef- 
fective mesons. Since it is well know since the early days 
of Skyrme theory that pure ^-forces are not sufficient to 
describe at the same time nuclear binding energies and 
radii, and since gradient terms in the Lagrangian can 
lead to certain difficulties in the relativistic formulation, 
historically the second method was the first to be used 
d, H^. Only recently relativistic point coupling mod- 
els with density dependent coupling constants have been 
employed successfully in nuclear physics [2X] . 

For simplicity we concentrate in this work on meson 
exchange models with non-linear meson couplings. Of 
course, the corresponding equations can be easily ex- 
tended to meson coupling models with density depen- 
dent vertices [28l . [29l . [soj or to relativistic point coupling 



models 

In covariant density functionals with meson exchange, 
the nucleons are described by Dirac spinors coupled by 
the exchange of mesons and by the electromagnetic field 
through an effective Lagrangian. The starting point for a 
phenomenological ansatz is therefore the Walecka model 
[23| . The mesons are classified by there quantum num- 
bers, spin, parity and isospin {I'^,T). In the isoscalar 
channel one has the scalar cr-meson {I^ ~ 0^,T = 0), 
and the vector w-meson {I^ = 1 ,T = 0), and in the 
isovector channel one considers only the vector p-meson 
{r = 1^,T = 1). The (5-meson {r = 0+,T = 1) is 
not included because, so far, there is not enough data in 
low energy nuclear structure physics to fix its parameters 
uniquely. In addition, the pion is not taken into account 
because, again for the sake of simplicity, we work only 
at the Hartree level, which forbids the appearance of the 
parity violating pion-field. The essential contributions 
of pionic degrees of freedom by two-pion exchange are 
taken care of in a phenomenological way by the tr-meson. 
Therefore, the starting point is an effective Lagrangian 
density of the form 

£ = Cn + £m + Cint- (1) 

Cn refers to the Lagrangian of the free nucleon 

Cn = ^ih^'d^ - m)^, (2) 

where m is the bare nucleon mass and ^ denotes the 
Dirac spinor. is the Lagrangian of the free meson 
fields and the electromagnetic field 

- \r,.R^'' + \mlp,(^ {f^.F^'" (3) 

with the corresponding masses m„, m^, rup, and the field 
tensors 

^fiu = d^uju — dijUjp^ 

Rfii^ = d^p„ - d„Pfj,, (4) 
Fpi, — d^Ajj — d^Ap, 

where arrows denote isovectors and boldface symbols 
vectors in 3-dimensional r-space. The interaction La- 
grangian Cint is given by minimal coupling terms 

Ant = -ga'4^T„a^l)-g^-ipT';^i^^'ilj-gp^T^tp^'ip-ei!V'^Af^'ip 

(5) 

with the vertices 

(6) 

where g^., g^, gp and e are the respective coupling con- 
stants for the (7, cj, p and photon fields. This yields to 

m 
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where the index m runs over the various meson and elec- 
tromagnetic fields, and also over the Lorentz index for 
vector mesons and isospin indices for mesons carrying 
isospin 



(r.,ri;,r;^,r^). (8) 



Already in the earliest applications of the RMF frame- 
work it was realized, however, that this simple linear in- 
teraction density functional did not provide a quantita- 
tive description of complex nuclear systems; an effective 
density dependence needs to be introduced. Historically, 
the first was the inclusion of non-linear self- interaction 
terms in the meson part of the Lagrangian in the form of 
a quartic a potential 



1 



2 2 



Uia) 



with 



rrf X 52 3 , 53 4 
U{<J) = yO-^ + —a , 



(9) 



(10) 



which includes the non-linear a self-interactions with 
two additional parameters (72 and 33. This particular 
form of the non-linear potential has become standard 
in applications of RMF functionals, although additional 
non-linear interaction terms, both in the isoscalar and 
isovector channels, have been considered over the years 

Two other approaches, of more recent development, 
can also be found in the literature, based on the intro- 
duction of the density dependence directly in the coupling 
constants [H, [2^ |30| and on the expansion of the meson 
propagators into zero-range couplings and gradient cor- 
rections terms f27j. For the sake of simplicity we will 
restrict the discussion in this work to non-linear density 
functionals, always taking the NL3 ^3B\ parameter set as 
the force of choice. Also, the explicit inclusion of the 
non-linear meson potential U{a) is generally avoided in 
order to keep the formulation as clean as possible. But, 
of course, it is included in all numerical calculations of 
results presented in the manuscript, and it shall be ex- 
plicitly mentioned at certain important points in the fol- 
lowing discussions. 

The Hamiltonian density can be derived from the La- 
grangian density of Eq. ([T|) as the (0,0) component of the 
energy-momentum tensor 



dqj 

leading the to the energy functional 
E[p,(l)] = / Hd^r. 



(11) 



(12) 



Following the Kohn-Sham approach 0, H^, one can ex- 
press the relativistic energy density i? as a functional of 
the relativistic single particle density matrix 



A 



(13) 



and the meson fields (f> = W^^iPil)- The sum over i in 
Eq. (|13p runs over the orbits in the Dirac sea {no-sea ap- 
proximation^ see below). Considering the 4-dimensional 
Dirac spinor ipi as a column vector and as a row vec- 
tor one concludes that p(r, r', t) is a 4x4 matrix in Dirac 
space. This leads to the standard relativistic energy den- 
sity functional 



Ermf[p, 4>\ = Tr[(-iaV + (3m)p\ + Tr[(/3r,„(/.„)/5] 
1 



(14) 



where summation over the different mesons is implied, 
and the trace operation involves summation over Dirac 
indices and an integral over the whole space. F™ de- 
scribes the structure of the meson-nucleon interaction. 
The upper sign in Eq. (fT4| holds for the scalar mesons 
and the lower sign for the vector mesons. At the mean 
field level the mesons are treated as classical fields. The 
nucleons, described by a Slater determinant |<i>) of single- 
particle wave functions, move independently in these 
classical meson fields. One can thus apply the classical 
time-dependent variational principle 



Jti 

which leads to the equations of motion 

idtp = [h[p,4>],p\ , 



J 0m = TTr [(3T„ip] , 



(15) 



(16) 
(17) 



where the single particle effective Dirac Hamiltonian h 
is the functional derivative of the energy with respect to 
the single particle density 

/i[/5,0] = ^^^ = (-*aV + /3m) + ^/?r„</.„. (18) 

m 

The time-dependent Dirac equation for the nucleons 
reads 



with the scalar S and vector potentials 
S{r,t) = g„a{r,t), 
Vi^ir, t) = g^ujf,{r, t) + gpfpf,{r, t) + eA^(r, t) 



(19) 



(20) 
I-T3 



and the time-dependent meson equations have the form 



[d^di, -f ml] a = -g^Ps 
[d-d,] = +ej^' 



(21) 
(22) 
(23) 
(24) 
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with the sources 



scalar — isoscalar ps — "01 "01 



vector — isoscalar 



vector — isovector 



electromagnetic 



A 



(25) 



(26) 



(27) 



A 1 



In order to describe the ground state properties of even- 
even nuclei, one has to look for stationary time-reversal 
invariant solutions of the equations of motion Eq. (jl6p 
and Eq. (|17p . The nucleon wave functions are then the 
eigenvectors of the stationary Dirac equation, 



-iaV + Va+ I3{m + S)] ipk = Ski'k 



(29) 



which yields the single particle energies £k as eigenvalues. 

The meson fields and the Coulomb potential obey the 
Helmholtz and Laplace equations 



9pPtv, 
epc, 



[-A] A° = 
with the following the source densities 



scalar — isoscalar 



vector — isoscalar 



i 

A 



(30) 
(31) 
(32) 
(33) 



(34) 
(35) 
(36) 
(37) 



Eq. (Ull) together with Eqs. jSOll-jSSl) pose a self- 
consistent problem which is readily solved by iteration. 
With the resulting density p and fields (p, the total en- 
ergy of the system can be calculated using Eq. (fH)) . Radii 
and other bulk properties of the nucleus can be derived 
as well. 

An important point of the present versions of covari- 
ant density functional theory is the no-sea approxima- 
tion, i.e. in the calculation of the sources for the meson 
equations ([30 |) -([33 |) only positive energy spinors are in- 
cluded in the summation. In a fully relativistic descrip- 
tion also the negative energy states from the Dirac sea 



vector — isovector ptv = E] V'J'ns'i/'ij 

i 

electromagnetic pc = E] V']^ ^ (1 ~ ■''3)V'i 



would have to be included. However, this would lead to 
divergent terms which have to be treated by a pro per 
renormalization procedure in nuclear matter [btI. l38l| or 
in finite nuclei [H, H^, |4l|, H^l • Numerical studies have 
shown that effects due to vacuum polarization can be as 
large as 20%-30%. Their inclusion requires a readjust- 
ment of the parameter set for the effective Lagrangian 
that leads to approximately the same results as if they 
were neglected [H, [H, This means that in a phe- 
nomenological theory based on the no-sea approximation, 
where the parameters are adjusted to experimental data, 
vacuum polarization is not neglected, it is just taken into 
account in the phcnomenological parameters in a global 
fashion and the no-sea approximation is in reality not an 
approximation. It is used in all successful applications of 
covariant DFT. This has, however, serious consequences 
for the calculation of excited states in the RPA [13 |. 

The vibrational response of the system can be stud- 
ied considering harmonic oscillations with small ampli- 
tude and with eigen-frequencies f2y around the station- 
ary ground p^'^\ In this case, the time-dependent density 
can be written as 



p{t) = p(0) -I- (,5p('')e-*"-* + h.c). 



(38) 



Imposing the condition that p is a projector at all times, 
the transition density matrices 6p^'^'^ have only matrix 
elements which connect occupied and unoccupied states 

[nl 



xl:I=Spl:l = {0\ala^W) 
0\al^a,\iy) 



mi rim 



(39) 



with respect to the stationary solution p^^\ In the non- 
relativistic case these are only ph- and hp- matrix ele- 
ments, i.e., the index i runs over all levels in the Fermi 
sea and the index m runs over all empty levels above the 
Fermi sea. 

In linear order, the equations of motion (|17p can be 
written down as the RPA equations in their standard 
matrix form 



A B 

-B* -A* 



(40) 



where the X^'^^ refers to the forward amplitude transi- 
tion density and Y^'^'> to the backward amplitude. The 
forward amplitude is thus associated with the creation 
and the backwards amplitude with the destruction of a 
p/i-pair. 

In the relativistic case the situation is more compli- 
cated. Because of the no-sea approximation in the RMF- 
model, the Dirac sea is empty. Therefore, we have to 
consider in relativistic RPA not only the ph- (and hp-) 
matrix elements of dp, but also the matrix elements Spah 
and Spha connecting states in the Dirac sea with those in 
the Fermi sea. The index i in the amplitudes and 
, and therefore in the RPA equation , runs again 
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over all the levels in the Fermi sea; however, the index 
m runs now over all the levels above the Fermi sea and 
over all the levels in the Dirac sea. This means we have 
not only to take particles in the Fermi sea and put them 
in the empty levels above the Fermi surface, but we have 
to consider also configurations where we form holes in 
the Fermi see and occupy empty levels in the Dirac sea. 
At a first glance this seems to be completely unphysical, 
because according to Dirac, the Dirac sea should be filled 
with particles. It turns out, however, that this is not the 
case. Considering the time-dependent RMF equations, 
the Dirac sea depends on time and the no-sea approx- 
imation should be realized at every point in time. In 
fact, solving these equations we consider only the time- 
evolution of the levels tpi{t) in the Fermi sea. The cor- 
responding time-dependent levels in the Dirac sea stay 
empty for all times [12]. When we describe this situation 
in the static basis, it is a mathematical consequence of 
the completeness of the basis, that one has to include 
also the a/i-configurations. If one neglects those config- 
urations, self-consistency is violated and one does not 
preserve the nice properties of RPA, such as current con- 
servation [455 and the separation of the Goldstone modes 
(spurious states) from the other physical solutions. 

Neglecting the a/i-configurations in the RPA equations 
also leads in specific cases to highly unphysical results, 
as for instance shifts in the energy of the GMR in ^"^Pb 
from the experimental value at 14 MeV down to 2-3 MeV 
[3. 

Taking into account also a/i-configurations renders the 
solution of the relativistic RPA equations much more 
complicated than in the non-relativistic case, because 
(i) the dimension of these equations increases consid- 
erably and (ii) the matrix A ± B is no longer positive 
definite and therefore the non-Hermitian matrix diago- 
nalization problem can not be transformed into a Her- 
mitian problem of half dimension, as discussed for in- 
stance in Ref. 11]. Only recently it has been shown that 
the relativistic RPA-equations can be reduced to a non- 
Hermitian diagonalization problem of half dimension [47| . 

For two different RPA excited states v and 1/' , the fol- 
lowing orthogonality relation holds 

^mi ^mi ^ mi ^ mi — "I'f ' J \^^J 

which can be used to normalize the eigen vectors 
(X^"^), F^"^)). Within the RPA-approximation the transi- 
tion matrix elements for a one body operator O between 
the excited state \v) and the ground state ]0) are given 

by 

(oiajz.) = Y.'^m^xi:] + 01,^^. (42) 

mi 

The RPA matrices A and B read 

Ami,nj — {Sm — £i)SmnSij + ^ji„, (43) 
Bmi,nj = Vmnij' (44) 



where the matrix elements VJ^^/j,/; are the second deriva- 
tives of the energy functional with respect to the single 
particle density 

ykf,'i = {knV^''\k'l)^T;f^ (45) 

OPk'kOpw 

and yP'^ is the effective interaction. As we have seen, 
the mean-field ground state is characterized by the sta- 
tionary density matrix p*-'^-' and by the meson fields , 
which, up to this point, have been treated as indepen- 
dent variables connected to the density by the equations 
of motion in Eq. (fT7|) . 

In order to describe small oscillations self-consistently 
it turns out to be useful to eliminate the meson degrees 
of freedom from the energy functional, such that the 
fermion equation of motion (approximated by the RPA 
equation (|43|) ) is closed, i.e., the residual interaction has 
to be expressed as a functional of the generalized density 
p only. This elimination of the meson degrees of freedom 
is possible only in the limit of small amplitudes, 

p = + dp, (46) 

where 6p and S(j) are small deviations from the ground 
state values p^*^^ and 0'-°-' . Substituting this expansion in 
the Klein-Gordon equations (|17p and retaining only the 
first order in Sp we find 

[d^d'' + mlJ] 6(l),n = T9mSpm, (47) 

with the local densities, the sources for the 
various meson fields are given by (5pm(r) = 
{Sps{r),Spy{r),Sp^t{r),Spc{r)) for m = {a,uj,p,A). 
Neglecting retardation effects (i.e. neglecting df) one 
finds for the linearized equations of motion for the 
mesons 

[-A + m^)] (50™ = TgmSpm- (48) 

This approximation is meaningful only at small energies, 
as compared to the meson masses, where the short range 
of the corresponding meson exchange forces guarantees 
that retardation effects can be neglected. A formal solu- 
tion for can be written as 

S(bmir) = tJ d\'g^G^{r,r')5p„,{r') (49) 

which allows us to decompose the residual interaction 
yph jj^ various meson exchange forces 

VP'' = V, + V^ + Vp + V-, (50) 

with 

V;,(l,2) - T<?^(/3r„)«G„(ri,r2)(/?r„)(2). (51) 

For linear meson couplings the propagator Gm obeys the 
Helmholtz equation 

{-^ + mlf)Gm{r,r')^6{r-r'), (52) 
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and has Yukawa form 



Gm(ri,r2) 



I g-"lmkl-r2| 

47r |ri-r2| 



(53) 



The vertices f3Tm reflect the difl^erent covariant structures 
of the fields as defined in Eq. Combining the spa- 
tial coordinates r and the Dirac index a = 1 ... 4 to the 
coordinate 1 = {ri,ai) we can express the relativistic 
two-body interactions in the following way 

• for the cr-exchange 

F.(l,2) = (54) 



• for the cj-exchange 



^2 p-m,„|j-i-7-2 

VL(l,2) = f^(l-aWa(2))r 



4tt 



• for the p-exchange 



\ri - r2\ 



(55) 



V,il, 2) = |£Lf(i)f(^)(l - aW«(2))^ (56) 

47r ri - r2 



• for the electromagnetic interaction 
K„,(l,2) = - ^ ^ (57) 



47r 2 



- ■^21 



In the case of non-linear meson couplings the Klein- 
Gordon equation l48l is replaced by 



[-A + ml]a + U'{<t) = -g^Ps 



(58) 



Considering small oscillations around the static solution 
cr(°\ it leads, instead of Eq. (gS]), to 



[-A + ml + W{r)\ Sa - -gjps, (59) 



with 



W{r) = U"{a^°\r)) (60) 

and the propagator Ga{r, r') obeys the equation 

[-A + ml + W{r)] Ga{r, r') = 5{r - r') (61) 

which cannot be solved analytically. More details, how 
to determine this propagator numerically are given in 
Appendix [B] 



concept of nuclear shape vibrations. If a system is de- 
formed, its spatial density is anisotropic, so it is possible 
to define its orientation as a whole, and this naturally 
leads to the presence of collective rotational modes. In 
1950, Rainwater [i^ observed that the experimentally 
measured large quadrupole moments of nuclei could be 
explained in terms of the deformed shell model, i.e., the 
extension of the spherical shell model to the case of a de- 
formed average potential. In a following paper [s^l , Age 
Bohr formulated the basis of the particle-rotor model, 
and introduced the concept of an intrinsic (body-fixed) 
nuclear system defined by means of shape deformations 
regarding nuclear shape and orientation as dynamical 
variables. The basic microscopic mechanism leading to 
the existence of nuclear deformations was proposed by 
A. Bohr [5l|, stating that the strong coupling of nuclear 
surface oscillations to the motion of individual nucleons 
is the reason for the observed static deformations in nu- 
clei. Nowadays, the deformation mechanism in nuclei is 
well understood [l2|: for sufficiently high level density 
in the vicinity of the Fermi surface, or for sufficiently 
strong residual interaction, the first 2+ excited state (a 
quadrupole surface phonon) is shifted down to zero en- 
ergy (it "freezes out"), effectively creating a condensate 
of quadrupole phonons and such giving rise to a static 
deformation of the mean-field ground state. 

In order to calculate excitations in deformed nuclei 
RPA theory outlined in the previous section can be used. 
It is important to remember, however, that these excita- 
tions are intrinsic in as much as they are relative to the 
local deformed ground-state. Nevertheless, the applica- 
tion of RPA for the calculation of intrinsic excitations of 
deformed nuclei is formally completely analogous to that 
for spherical nuclei. The only difference is that one has 
now single particle orbitals violating rotational symme- 
try, i.e. having no good angular momentum. For this 
reason it is not possible to apply group theory and to 
reduce the dimension of the RPA matrix by angular mo- 
mentum coupling techniques. Only in the case of axial 
symmetry, reductions are possible which using the good 
quantum number K, the projection of the total angular 
momentum on the symmetry axis. 

The introduction of a deformed intrinsic state in DFT 
is straightforward. Let us suppose that there exist a sym- 
metry operator S such that the energy density functional 
is invariant under the symmetry transformations e*"*^, 
i.e. for a transformed density pa 



we have 



E[p] 



(62) 



(63) 



III. RMF+RRPA IN DEFORMED NUCLEI 

The fact that nuclei can be deformed was already em- 
phasized by Niels Bohr in his classic paper on the nu- 
clear liquid-drop model [481] , where he introduced the 



Examples of such a symmetry in even-even nuclear sys- 
tems would be rotational and translational symmetries 
and the third component of isospin (i.e. the charge). If 
the density has the same symmetry, pa = p-, we can re- 
strict the set of variational densities to those with this 
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symmetry. However, such a symmetric solution is not 
necessarily at the minimum in the energy surface defined 
by E[p], that is, the best solution. Because of the non- 
linearity of the variational Eq. (|29p it is possible that the 
solution breaks the symmetry spontaneously, i.e., the en- 
ergy density is invariant under 5-transformations but the 
density is not pa ^ p- 

Rotations are one of such continuous symmetry trans- 
formations. Nuclei with semi-closed shells have pair- 
ing correlations, and the solution with lowest energy of 
the variational Hartree-Bogoliubov (HE) or Hartree-BCS 
(HBCS) equations have a spherical intrinsic density dis- 
tribution. One can always write their ground state wave 
function as a rotationally invariant product state of HB 
or BCS type . On the other hand, most nuclei through- 
out the periodic table have open shells for both types 
of particles, and thus due to the strongly attractive se- 
niority breaking proton-neutron interaction their respec- 
tive intrinsic single particle densities are usually not in- 
variant under rotations. Nevertheless most of the nuclei 
have minima with axial symmetric density distributions. 
There are only few cases with pronounced triaxial defor- 
mations. 

The present investigation is restricted to nuclei that 
can be adequately described by a variational wave func- 
tion with axial symmetry, and so the projection of the 
angular momentum f2 on the symmetry axis is a con- 
served quantity. It is therefore convenient to use cylin- 
drical coordinates 



r = (r cos if, r sin ^p^ z)^ 



(64) 



were, as usual, the symmetry axis is labeled as the z- 
axis. Note, that r is the distance from the symmetry 
axis, not the distance from the origin. For reasons of 
simplicity we avoid in this manuscript, as far as possible, 
the notation r^. The single-particle Dirac spinors 'lpk^ 
solution of Eq. , are then characterized by the angular 
momentum projection $1, the parity tt and the isospin 
projection i, 



/ /+(r,z)e^("»-V2)0 \ 
Vifffc(r,z)e*(f^'+i/2)<A j 



(65) 



For even-even nuclei, for each solution tl)^ with positive 
il/c there exists a time-reversed one with the same energy 
that will be denoted by a bar 



(66) 



The time reversal operator has the usual form ij^j^K, 
where K is the complex conjugation. 



A. Configuration space for the RPA equation 



be formed using the single-particle spinor solutions of 
the static problem. Since the total angular momentum is 
no longer a good quantum number, we cannot take ad- 
vantage of angular momentum techniques when forming 
these pairs. Only axially symmetry and parity is left. 
The full RPA-matrix can be thus reduced to blocks with 
good quantum numbers K and tt. In particular, this 
means that the RPA matrix elements V^^j„ must obey 
the following selection rules 



Thus we can define RPA phonon operators 



'v,K'" / y ^^mi "-m^i ^ mi "'i^m 



(67) 
(68) 



(69) 



as linear combination of pairs with good angular momen- 
tum projection K and parity tt. This mean that the sum 
runs only over pairs mi such that the conditions (|67p and 
(l68t are satisfied and the different excitation modes 



W,Kn=Q+K-\0) (70) 
can be labeled by the quantum numbers 

/C'^ = 0±,1±,2=^,--- (71) 

where 

- (f7„, - f],)^"'""'^ (72) 

One has to be careful handling time reversal symmetry 
in the case of coupling to K = 0, where for each pair of 
the form of Eq. (|72|) there exists the time reversed one 



= i-n^ + 17,)("""') = (r!„ - r!,)^"""*) (73) 



with the same energy that also satisfies Eq. ([67]) , and has 
to be considered explicitly when calculating the matrix 
elements. 



B. Evaluation of tlie RPA matrix elements 

As we have seen in Eq. (P5|) . the matrix elements of 
the residual interaction can be derived from the energy 
functional as the second derivative with respect to the 
density. In the case of meson exchange models this in- 
teraction in Eqs. ([50]) and ([?T|) . The index m runs over 
the various mesons, but, for vector mesons, also over the 
Minkowski index p. For simplicity we therefore use in 
the following for the vertices only the Minkowski index 
p instead of m. In this case the interaction has the form 

l4(l,2) =Tg^/3Wr^«G,„(ri,r2)/3(2)r(f) (74) 



The rows and columns of the RPA matrix in Eq. (HD)) 
are labeled by all the possible ph- and a/i-pairs that can 



where the propagator G{ri,r2) has Yukawa form for 
mesons with linear couplings, and has to be evaluated 
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numerically in the other cases. We first concentrate on 
mesons with linear couplings. In this case the propagator 
Gm{fi — ^2) depends only on ri — r2 and can be written 
in Fourier space as 



)e- 



with the meson propagator 



Am(q) 



q2 + rn^ 



The interaction ([7H) has the form 



with 



(75) 



(76) 



XI, 2) = t/ (^0^(9- l)^MQl{<l, 2) (77) 



(78) 



mi.nj — (^m ^i)^mn^ij~^ 



the RPA matrix 1140 
(2^ 



(27r)3 



m|g'^(q)|z)A„(q)(n|g^(q)br (81) 
;m|Q'^(q)|i)A„(q)(J|QM('7)l")* (82) 



where \n) = T\n) is the time- reversed state to |n). Us- 
ing the symmetry properties of the operators Q'^{q) one 
obtains 



Q\QM)\^) = {~)Hn\QM)\3) 



(83) 



where S is the spin of the exchanged meson, i.e. 5 = 
for scalar mesons and the time-like part of vector mesons 
and S = 1 for the spatial part of the vector mesons. 



C. Matrix elements in the intrinsic and in the 
laboratory frame 



For each q, Q^{q) is a one-body operator in r-space and 
in the 4-dimensional Dirac-space defined by the combined 
index 1 = (ri, ai). This definition of the operator is 
flexible enough to allow also applications of the meson ex- 
change model with density dependent coupling constants 
gmix) = gm{p{,i'))- In the present investigation, however, 
we do not follow this avenue. Considering the q-integral 
as a sum over discrete values in g-space, the interaction 
(177)) is a sum of separable terms. The corresponding two- 
body matrix elements can thus be expressed by the one- 
body matrix elements of the operators Q^{q)- Using this 
form we can evaluate the two-body matrix elements 

{ki'\vi^\k'i) ^ I ^(k\Q^iq)\k')AM{i\Q,iq)\ir 

(79) 

with the single particle matrix elements 

(fcig^(q)lfc') = J d\ graMr)T^{r)e"i''i^k'{r). (80) 

For the case with axial symmetry, the evaluation of these 
matrix elements is best accomplished in cylindrical coor- 
dinates . In this case one finds that the integrals over 
the azimuth angles in coordinate and momentum space 
can be evaluated analytically. This leads to the selection 
rule rife — ilfe' = ri; — f2;' (details are given in Appendix 

For non-linear meson couplings the propagator 
Gm{r,r') depends on both coordinates and therefore we 
find in Fourier space the matrix Am(<7, q'), which is cal- 
culated numerically by matrix inversion. This leads to 
a four-fold integral in momentum space for the evalua- 
tion of the two-body matrix elements (|79p (for details see 
Appendix [BJ . 

Summarizing this section, one finds for the elements of 



So far we have solved the relativistic RPA equations in 
the intrinsic frame. Neither the basis states, the p/i-states 
based on a deformed ground state, nor the eigenstates 
\i^,K'^) of the RPA equations in Eq. (|7D|) are eigenfunc- 
tions of the angular momentum operators and Jz in 
the laboratory frame. We therefore call the states {ly, K^) 
wave functions in the intrinsic frame. In fact, they have 
little in common with the wave functions in the labora- 
tory frame, which have to be eigenstates of the angular 
momentum operators and Jz- In order to calculate 
matrix elements which can be compared with experimen- 
tal data we therefore have to project onto good angular 
momentum, i. e. on eigen spaces of these operators 
and Jz in the many-body Hilbert space. 

Using the projection operators defined in Ref. [ll[ we 
obtain the wave functions 

\iy,K,IM) ^PIikW^K^) (84) 

where V{'}j^{fl) are the Wigner functions [52], the irre- 
ducible representations of the group 0(3) of rotations in 
3-dimensional space. They depend on the Euler angles 17, 
and R{fl) is an operator which rotates the intrinsic wave 
function K'^) by the Euler angles fl. The evaluation of 
matrix elements in the many-body Hilbert space using a 
projected wave functions is a rather complicated task. In 
involves in particular the calculations of the overlap in- 
tegrals {iy,K''\R{n)OR{n')y,K"'). it has been found, 
that for well deformed intrinsic wave functions these over- 
lap integrals are sharply peaked at — ft' . Replacing 
this sharply peaked functions by Gaussians with a rather 
small width and in the extreme limit of strong deforma- 
tion by S{Q — n') one obtains the so-called needle approx- 
imation [li| . The overlap functions are sharply peaked, 



9 



in particular for systems with many particles, and there- 
fore the needle approximation is not only valid in cases of 
strong deformations in the geometrical sense, but in gen- 
eral for heavy systems with normal deformations, and in 
the classical limit even for a spherical shape with a well 
defined orientation. Moreover, it can be shown that the 
results obtained with this approximate projection (the 
needle approximation), are e quiv alent to the results of 
the particle plus rotor model [51[ where the orientation 
n of the intrinsic frame is used as a dynamical variable 
(for details see appendix [E)) . 

The evaluation of the matrix elements in the labora- 
tory system reduces to the calculation of products of spe- 
cific intrinsic matrix elements and geometrical factors. 
This leads to the following expression for the reduced 
matrix elements 



{IfKf\\dx\\hK,)^i2h + l){2If + l) 



(85) 



I^ ^ If 

K, /i Kf 



+ (-1) 



{Kf\dx^\K^, 



where {Kf\0\^\Ki) is the intrinsic matrix element of the 
multipole operator O^/i which is easily calculated with 
the help of Eq. (gH). 

In the following we therefore have to distinguish matrix 
elements and transition densities in the intrinsic frame 
calculated directly with the solutions of the RPA equa- 
tion and matrix elements and transition densities in the 
laboratory system, which are obtained after angular mo- 
mentum projection in the needle approximation in Eq. 

(ii. 



IV. STRENGTH FUNCTIONS AND SUM 
RULES 

Experimental nuclear spectra show in the continuum 
excitations as resonances with finite width. Since the di- 
agonalization of the RPA equations is done in a discrete 
basis, we obtain discrete eigenstates \ v). Using Eqs. (|42l) 
and we can calculate for each of them the reduced 
transition matrix elements for specific multipole opera- 
tors, as for instance the reduced transition probabilities 
{BEI- and i?Af/-values) for electric and magnetic tran- 
sitions 



['^JWQikWO) (86) 



B{MI,0^ I,K,uj,) = {i^JWMikWO) 



(87) 



It is well known that the width cannot be described well 
within the RPA approach discussed here. On one side we 
work in a discrete basis, and therefore the continuum is 
not treated properly and the escape width is not taken 
into account, and on the other side RPA itself is a lin- 
ear approximation. It does not contain the coupling to 



2p2h- and more complicated configurations and therefore 
does not allow a proper treatment of the decay width. 
Higher order correlations, for exam ple the coupling to 
low- lying collective phonons [H, Is^. ISSj , have to be in- 
cluded for this purpose. It is, however, also known from 
spherical RPA calculations, that this method is able to 
describe rather well the position of the resonances and 
the strength of the transitions for given multipole oper- 
ators, i.e. the percentage of the sum-rule exhausted by 
a specific resonance. To overcome the problem of the 
width, we adopt a phenomenological concept and aver- 
age the discrete RPA strength distribution obtained from 
the solution of the RPA equations in a discrete basis with 
a Lorentzian function of a given width P. For the electric 
response we have: 



R{E) B{EI,0^ I,K,uj,)- 



r/2 



(r/2)2 



t: (E — LUu 



and a corresponding expression holds for the magnetic 
response. This results in continuous strength functions 
which can be compared with experimental spectra and 
sum rules. The knowledge of the sum rules is of special 
interest, since they represent a useful test of the models 
describing the collective excitations [ll|. For example, 
the energy weighted sum rule (EWSR) for a transition 
operator O can be represented as a double commutator 



^1 



(89) 



If one assumes a non-relativistic Hamiltonian, a local op- 
erator O and a local two-body interaction, only contri- 
butions from the kinetic energy contribute and one can 
evaluate this sumrule in a model independent way: 



2m 



An 



(90) 



These classical values for the sum rules are only approx- 
imate estimates. In practical calculations they may be 
enlarged by an enhancement factor due to the velocity 
dependence and due to exchange terms of the nucleon- 
nucleon interaction. It can be shown, that many of these 
sum rules apply also in RPA approximation. In this 
manuscript We evaluate the EWSR in the interval be- 
low 30 MeV excitation energy as 



(91) 



In particular, the Lorentzian function in Eq. (1881) is nor- 
malized in such a way as to give the same EWSR as 
calculated with the discrete response 

Si = Y^ LU^B{EI, 0^ I,K,LU^)= j E R{E)dE (92) 



Sum rules also offer the possibility of a consistent defi- 
nition of the excitation energies of giant resonances, via 
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the energy moments of the discrete transition strength 
distribution 



(93) 



In the case k = 1 this equation defines the energy 
weighted sum rule of Eq. ([91]) . If the strength distri- 
bution of a particular excitation mode has a well pro- 
nounced and symmetric resonance shape, its energy is 
well described by the centroid energy 



E = 



Alternatively, mean energies are defined as 



mk-2 



(94) 



(95) 



where the difference between the values Ei and E^ can 
be used as an indication of how much the strength distri- 
bution corresponding to an excitation mode is actually 
fragmented. If the multipole response is characterized by 
a single dominant peak, the two moments are equal, i.e. 
El = E3. In the relativistic approach, due to the no-sea 
approximation, the sum in Eq. (|93p runs not only over the 
positive excitation energies, but also includes transitions 
to the empty states in the Dirac sea which contribute 
with negative terms to the sum. As it was pointed out 
in [5I, [53, HI] , for the EWSR the double commutator of 
Eq. ([5^ should vanish, and it is another good check for 
the numerical implementation. 



V. TRANSITION DENSITIES 

In order to have a intuitive picture of the nuclear exci- 
tations we investigate in the following the time evolution 
of the baryon density. Let us consider the baryon four- 
current operator in coordinate space 



(96) 



with single-particle matrix elements in the Dirac basis 

Jkk'ir)=Mr)rA'{r), (97) 
which can be written as 



J' 



ak' 



(98) 



fcfc' 



In order to calculate its time evolution within the RPA 
approximation, and for a particular excitation mode v, 
we use Eq. and find 

S.fir) = Y.(^L(r)X^ + jLirTY^:^) (99) 



Thus, the total time dependent baryon four-current for 
a given excitation mode v with energy lo^, is 

rir, t) = j^ir) + <5j^(r)e-''^'^* + <5j^(r)*e'-'^* (100) 

In particular, the baryon density p{r, t) = j'^{r, t) can be 
written as 



p(r, t) = po(r-) + 5p(r)e-*'^'^* + 5p(r)*e'' 



(101) 



Throughout the rest of this paper, all types of intrin- 
sic transition densities refer to the baryon intrinsic tran- 
sition density in coordinate space, 5p{r), as defined by 
the zero's component of Eq. As discussed in the 

previous section, we also have to distinguish the transi- 
tion densities in the intrinsic system from those in the 
laboratory frame. 

In a classical system the transition density would de- 
scribe the actual movement of particles. In the quantum 
mechanical description one considers a time-dependent 
wave packet and decomposes it into the contributions of 
the different excited eigenstates of the system. The tran- 
sition density is an off-diagonal matrix element between 
the stationary ground state |0) and the excited eigen state 
I v) and it is regarded as a measure of the contribution of 
the eigen state \v) of the system to the evolution of the 
time-dependent wave packet. To what extent an state 
\v) can be interpreted in the classical sense depends on 
percentage of the sum rule exhausted by the transition 
strength of this excitation mode Therefore the tran- 
sition densities provide an intuitive understanding of the 
nature of the excitation modes. It can be used as well for 
the calculation of transition probabilities as for a quali- 
tative understanding of these modes. 

In an axial symmetric system the transition density 
in (IMl) can be written as 



5p{r) ^dp{r^,z)e-'^'^ 



(102) 



where K is the angular momentum projection of the ex- 
citation mode under study. Note, that we distinguish in 
this section the coordinates r± (the distance from the 
symmetry axis) and r (the distance from the origin). 
Substituting this last expression in Eq. (jlOip we arrive 
at 



p{r±,^, z, t) = po(r±, z) + [jp(ri, z)e-*(^'^+'^''*) + h.c 

(103) 

The two dimensional quantities Sp{rj_,z) will be plot- 
ted when discussing intrinsic transition densities, and 
no further reference to the phase expressed in the ex- 
ponentials will be made. To interpret these plots, it 
is useful to keep in mind that 6p{rj_,z) has to be con- 
sidered together with Eq. (|103p in order to obtain the 
full three dimensional geometrical picture. However, to 
be able to compare with experimental transition densi- 
ties measured in the laboratory frame, we project the 
two dimensional intrinsic transition densities 5p{r±,z) 
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on to good angular momentum. This is done by expand- 
ing the current operator in Eq. (|98p in spherical coordi- 
nates r — {r sin 9 cos (p,r sin 9 sin ip,r cos 9) using the set 
of spherical harmonics Ylm(^<p) as a basis 

fir)= J2 jlkk'ir)YLK{0v) alak' (104) 

LK,kk' 

where 

JL,fcfc'(0 = / dcos9d^ji:,,{r)Ylj,{9^). (105) 

The projected transition density reads 

6pir)=SpLir)YLK{9^) (106) 
with the radial projected transition density 

SpLir)^ J dcos9dipSp{r^,z)YlK{9if ) (107) 

In the following figures we present the quantity r^SpL- 
Because of the approximations in the derivation of 
Eq. l|55|) . this last equation only holds approximately. 
Nevertheless, we will see that the results for well de- 
formed nuclei are excellent. For example, the transition 
density patterns for the Giant Dipole Resonances (GDR) 
and the Pygmy Dipole Resonances (PDR), are in reason- 
able agreement with those found experimentally and in 
other theoretical RPA studies in spherical symmetry. 

VI. SPURIOUS MODES AND NUMERICAL 
IMPLEMENTATION 

If the generator for a symmetry operation of the full 
two-body Hamiltonian, which is represented by a one- 
body operator, does not commute with the ground state 
density, there exists a Goldstone mode, a so-called spuri- 
ous solution, of the RPA equations with zero excitation 
energy associated with this symmetry. These solutions 
are not really spurious, but they correspond to a collec- 
tive motion without a restoring force ^5£] and therefore 
they do not correspond to oscillations with small ampli- 
tude. Examples are translations or rotations. In princi- 
ple this should not be significant, since we are concerned 
only with the intrinsic structure of the nucleus. Thouless 
found that for the exact solution of the self-consistent 
RPA equations these spurious modes are orthogonal to 
all the other modes. They do not mix with them and can 
be separated [60l] . 

In practical applications, however, in many cases the 
spurious solutions are not completely orthogonal to the 
physical states by various reasons. One should be able to 
distinguish them from the true vibrational response of the 
nucleus, as experience shows that this mixture can lead 
to seriously overestimations in the strength distributions. 

In normal calculations, because of numerical inaccu- 
racies, truncation of the p/i-space and inconsistencies 
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FIG. 1: (Color online) Dependence of the — l""" rota- 
tional spurious mode on the coordinate and momentum mesh 
size for the non-linear model NL3 parametrization. For a co- 
ordinate and momentum mesh size of (41x41) and (31x31), 
respectively, the accuracy limit of the diagonalization proce- 
dure is achieved. The logarithmic scale in the z-axis is used 
to enhance the readability of the graph. The lowest z value 
corresponds to a value of around 0.05 MeV. 

among the ground state and RPA equations, the spuri- 
ous states are often located at energies somewhat higher 
than zero and often cause a mixing with physical states. 
There are several approaches to overcome this problem. 
Some authors adjust a free parameter of the residual in- 
teraction until the energy of the spurious mode goes to 
zero. Another method is to remove a posteriori the spu- 
rious components from the physical states by projection. 
This is possible, because the wave functions of the spu- 
rious modes are given by the matrix elements of the cor- 
responding generators |ll| . 

In this investigation a fully self-consistent implemen- 
tation of the RPA is used, and thus as long as numer- 
ical inaccuracies are kept to a minimum, the spurious 
modes decouple without further complications. Because 
the block-wise structure of the RPA matrix, they are ex- 
pected to be present only for specific quantum numbers 
when specific symmetry constraints are met; since we are 
restricting to axial symmetry, their expected appearance 
can be summarized as 

• A rotational spurious mode for the = 1+ chan- 
nel associated with rotations of the nucleus as a 
whole around an axis perpendicular of the symme- 
try axis in z-direction. Its generator is the angular 
momentum operator = + iJy [6l|. 

• A translational spurious mode for the = 0^, 1^ 
channels associated with the translation of the nu- 
cleus as a whole. Its generators are the linear mo- 
mentum operators Pz and P-^- . 

In fact, the position of the spurious modes provides a 
very accurate test of the actual numerical implementa- 
tion of the RMF+RRPA framework. Thus, it is impor- 
tant to study their evolution with the approximations 
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FIG. 2: (Color online) Spurious — 1"*" rotational mode 
dependence on the maximum interaction momentum while 
keeping the number of mesh points constant. Good numerical 
results for a momentum mesh size of (31x31) can be achieved 
with a maximum momentum in the interval 5 < q,nax < 9. 



performed. In the present status of the implementation 
seven parameters control the numerical accuracy, and can 
be categorized in two groups. The first group specifies the 
precision of the numerical integrations. In this category 
are included the number of coordinate and momentum 
lattice points and the upper boundary of the momentum 
integrals. The second group deals with the size of the 
configuration space, and includes the energy cutoffs for 
ph- and a/i-pairs. 

Since it is unfeasible to study this seven-dimensional 
surface in detail, when studying the dependence of the 
spurious modes on one, or a set of, parameters, those not 
under scrutiny were fixed to best possible values the hard- 
ware supports. This means, in particular, that the full ph 
configuration space is taken if not otherwise stated, and 
that the maximum momentum is fixed to gmax — 8 fm"-'^, 
well above the Fermi momentum of the nucleus. 

In Fig. [T] the position of the rotational K"^ = 1+ spuri- 
ous mode in ^'^Ne is plotted against the number of points 
in the coordinate and momentum lattices. For a rela- 
tively low number of points a plateau is reached where 
further improvement of the accuracy cannot be achieved. 
The optimal number of evaluation points for the quadra- 
tures is therefore around 41x41, which allows for very 
precise calculations. Furthermore, additional tests show 
that the overall precision in the determination of the 
energy of excited states of the code is capped out at 
0.05 MeV, which is surprisingly good. 

Figure[2]depicts the position of the rotational K"^ = 1+ 
spurious mode for ^"Ne against the maximum momentum 
of the expansion used in the integral for the evaluation 
of the single particle matrix elements in Eq. ([50)) . The 
flat region between 5 and 9 fm~^ hints that a maximum 
momentum of gmax = 5 fm~^ provides enough precision 
for the proper decoupling of the spurious mode. The 
increase observed in the position of the spurious mode 
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FIG. 3: (Color online) On the upper plot we show the de- 
pendence of the — 1^ rotational spurious mode on the 
size of the ph- and ah-space size. The lower plot shows the 
same dependence for the K'^ = 1"*" and = 0~ translational 
spurious modes. 



for maximum momentum values larger than 9 fm~^ is an 
artifact due to the number of points for the momentum 
lattice being fixed at (31x31). 

Figure [3] shows the dependence of the rotational = 
1"^ and translational = 0~, 1~ spurious modes on the 
configuration space size for ^°Ne, as calculated with the 
NL3 parameter set. In the translational case two curves 
are plotted, one for the K'^ = 0~ mode and one for the 
K"^ = 1~ mode. It is interesting to note that, even if 
the spurious mode can be brought very close to zero, it 
requires the inclusion of almost all the possible p/i-pairs in 
the configuration space. In this specific case, i.e. in ^''Ne, 
that amounts to the inclusion of roughly five thousand 
pairs. Several test indicate that the situation improves 
greatly in heavier nuclei, where usually 5% percent of 
all possible p/i-pairs are enough to decouple the spurious 
modes at energies around 0.5 MeV. 

Since for the solution of the ground-state the equations 
of motion are expanded in a harmonic oscillator basis, the 
configuration space where the RPA is solved does not 
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FIG. 4: (Color online) Dependence of the A"" = 1" and K'^ = 
0~ translational spurious modes on the configuration space 
size, dictated by the number of oscillator shells used in the 
ground-state calculation. 
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FIG. 5: (Color online) Response in Ne to the operator J±i, 
generator of rotations around a perpendicular of the symme- 
try axis. Almost 100% of the strength is exhausted by the 
spurious mode (situated at 0.08 MeV), with minimal admix- 
ture to the physical states. 



spawn the whole Hilbert space, even if all possible ph- 
pairs are taken. The quality of this expansion depends 
on the number of major oscillator shells used, and the 
results obtained at the RPA level will be influenced by 
this approximation. In particular, the proper decoupling 
of the translational spurious mode is very sensitive to the 
number of oscillator shells employed. In Fig.[4]the trans- 
lational spurious mode is plotted versus the number of 
oscillator shells used in the ground state. Already with 
the inclusion of 16 major shells is enough to achieve a 
precision in the spurious mode of around 0.1 MeV. In all 
practical cases presented in the this study the number of 
oscillator shells was chosen between 12 and 16, depend- 
ing on the desired final precision and the availability of 
computer resources. 

However, more important that the actual position of 



the spurious modes is their admixture to the real physical 
states. For the same reasons the spurious mode does not 
appear at exactly zero energy, the physical states are not 
completely orthogonal to it, producing unreal results and 
very often overestimated strength. 

Moreover, there is one important property of the spu- 
rious modes that can also be used to measure the extent 
of their admixture with the rest of the RPA states. They 
are not normalizable in the sense of Eq. (|4T|l because 



(108) 



However, in all numerical implementations the relation 
(|108p is only approximately fulfilled because the spurious 
modes do not decouple exactly. How good the decoupling 
of the spurious modes is can be measured comparing the 
relative norms of the different eigen-modes. For an ap- 
proximate spurious mode labeled as (sp) it should hold 



that X 
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while for any other RPA mode by initial assumption, 
it holds that Xmi ^ Ymi, and thus 
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Our tests indicate that when the spurious modes are lo- 
cated below 0.5 MeV, the value of S in Eq. (|109p is at least 
three orders of magnitude smaller than the ones belong- 
ing to normal RPA modes. This is a very good indicative 
of the proper decoupling of the Goldstone modes. 

As an example of the low admixture of spurious com- 
ponents with the physical states. Figure [5] shows the re- 
sponse to the generator of the rotational spurious mode, 
the operator J±i, which represents rotations around an 
axis perpendicular to the symmetry axis. More than 
99.99% of the strength is exhausted by the spurious 
mode, which is located below 0.1 MeV. Similar results are 
obtained for the translational spurious modes in ■^"Ne. 

In general, it was observed that, if the position of the 
spurious mode is below 1 MeV, the strength function of 
the rest of the spectrum is mostly unaffected. The spec- 
trum in the low energy region, below 5 MeV, is, however, 
more sensitive to admixtures of the spurious modes; as a 
rule of thumb, the confidence limit in the position of the 
spurious mode, for a proper decoupling, has been consis- 
tently found around 0.5 MeV. 

There is still another test that can be devised in or- 
der to check the consistency of the whole framework, 
namely the conservation of spherical symmetry. Even 
though all the formulas are particularized to the case of 
axial symmetry, the interaction is rotationally invariant, 
so they should still be valid when a spherical ground- 
state is taken as basis for the RPA configuration space, 
i.e., they should preserve spherical symmetry. 
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running times and thus in the possibihty of detailed anal- 
ysis. With the optimal number of oscillator shells for a 
ground state calculation with full precision, the number 
of pairs never exceeds five thousand. Furthermore, be- 
cause the number of protons and neutrons is the same, 
switching off the electromagnetic interaction should give 
identical results for both protons and neutrons. Using 
this technique, very detailed checks can be carried out 
on the isospin part of the interaction, and its consistency 
can be further established. All these reasons make ^°Ne 
the ideal theoretical playground where to introduce the 
concepts that can later be used in the study of more com- 
plex systems. In the this section we shall present two 
sample applications for the well deformed nucleus ^°Ne. 



FIG. 6: (Color online) K'" = 0" and K'^ = 1" response to 
the El transition operator for the spherical nucleus ^"^O. 

In Figureiniis plotted the El excitation strength for the 
spherical nucleus ^^O. Since the El operator is a rank-one 
tensor, it has three possible angular momentum projec- 
tions, K — that have to be calculated separately. 
The response for the modes with K — —1 and K = 1 are 
identical and correspond to vibrations perpendicular to 
the symmetry axis, i.e. one can calculate only one of 
them and double its contribution. The K = Q mode cor- 
responds to vibrations along the symmetry axis. If the 
nucleus is prolate, like ^°Ne, the response for in the K — Q 
mode should lie at lower energies than the K — 1 mode, 
as the potential is flatter in the direction of the symmetry 
axis. However, if the nucleus is spherical, like ^^O, there 
is no distinction between the K = and K — 1 modes, 
and their corresponding excitation strength should lie at 
exactly the same energies. From Figure [5] one can attest 
that the procedure for the solution of the RPA equation 
in axial symmetry indeed preserves rotational symmetry 
with a good degree of accuracy. 

The study presented concerning the decoupling of the 
spurious modes and the preservation of spherical sym- 
metry shows that the numerical implementation solves 
the equations posed by the self-consistent RMF+RPA 
framework in axial symmetry. We have also ascertained 
that a high degree of accuracy can be achieved in real 
calculations, as well as validated the good reproduction 
of formal and mathematical aspects of the RPA theory. 

VII. APPLICATIONS IN ^'NE 

As a first application of the RMF-f-RRPA framework 
we have undertaken a model study of the magnetic and 
the electric dipole response in ^°Ne. This nucleus offers 
several advantages. Its ground state is well deformed 
and exhibits a prolate shape in the RMF model, with 
a quadrupole deformation parameter [3 ~ 0.5. Another 
advantage is the reduced number of nucleons to be taken 
into account in the calculations, which translates in fast 



A. The magnetic dipole (Ml) response 

We first consider the magnetic dipole response. The 
discovery of low-lying Ml excitations, known as scissors 
mode, was made by Richter and collaborators in ^^^Gd 
in Darmstadt through a high-resolution inelastic electron 
scattering experiment [13] ■ The search for such a mode 
was stimulated by the theoretical prediction of a collec- 
tive mode, where the deformed proton distribution oscil- 
lates in a rotational motion against the deformed neu- 
tron distribution [63^, 'g^, [6^ Ha [s^ • The name "scissors 
mode" was indeed suggested by such a geometrical pic- 
ture. An excitation of similar nature was also predicted 
by group-theoretical models [gI, [gI, [tO] ■ 

The mode has been detected in most of the deformed 
nuclei ranging from the /p-shell to rare-earth and ac- 
tinide regions. The mode has been well characterized, 
and it is established that it is fragmented over sev- 
eral closely packed Ml excitations. For reviews to this 
mode and for recent semiclassical investigations see Refs. 

[lll[2i[2l,[2l,[2i,[2l 

A byproduct of the systematic study of the scissors 
mode was the discovery of spin excitations. Inelastic pro- 
ton scattering experiments on ^^^Sm and other deformed 
nuclei found a sizable and strongly fragmented Ml spin 
strength distributed over an energy range of 4-12 MeV 
[T^ . The experimental discovery stimulated theoretical 
investigations in the RPA approximation [T3]. 

Figure[7|shows the response to the Ml magnetic dipole 
operator in the nucleus ^°Ne . The shaded region corre- 
sponds to the full Ml response; the blue dashed line is the 
response to orbital part of the Ml operator and the red 
dotted line refers to the spin part. The calculations were 
performed with the maximum precision allowed by the 
current implementation of the computer code. The num- 
ber of pairs is around five thousand. Optimal numerical 
parameters were chosen in order to minimize the error. 
The rotational spurious mode is well separated, situated 
around 0.1 MeV, i.e. no admixture with the vibrational 
response is observed. 

Only one prominent peak is found around 5.7 MeV. 
Regrettably, there is no experimental data available for 
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FIG. 7: (Color online) Ml Excitation strength for ^°Ne, using 
the NL3 parameter set. A very well developed peak can be 
seen around 5.8 MeV. Its structure is composed mostly of spin 
flip transitions. 



the magnetic response in this nucleus. Theoretical stud- 
ies using large scale shell model calculations [ill predict a 
low lying orbital mode around 11 MeV for ^°Ne, in strong 
disagreement with our results. However, other calcula- 
tions [t^ performed in ^^Ne, with the same shell model 
interaction, exhibit two dominant low lying peaks around 
5-6 MeV. The orbital contribution to the total response 
is below 25%, which is in better agreement with results 
found within our RMF+RRPA calculations, where frag- 
mented strength with similar characteristics is found in 
the same energy region. 

Regarding the contributions from the orbital and spin 
components of the Ml operator to the total response, it 
can be observed in Figure [7] that there are two differen- 
tiated energy regions. Around the main excitation peak 
at 6 MeV there is an enhancement of the response due to 
the additive interference of the orbital and spin contribu- 
tions. In contrast, in the energy region above 6.5 MeV we 
observe the opposite, destructive interference and both 
contributions cancel. This feature of the Ml strength 
distribution has been also found in other studies [l3] and 
is much more evident in the case of heavier nuclei. 

From this figure we can also recognize that the main 
contribution to the total response comes from spin exci- 
tations. The supposed orbital character of the low lying 
spectra in the Ml transitions is eclipsed by the prepon- 
derance of spin flip strength, three times larger than the 
orbital response. Again, this is in disagreement with the 
cited shell model calculation [ts'], which in ^'^Ne predicts 
a much bigger orbital contribution to the total strength. 
However, low lying collective transitions in such a light 
nucleus as ^°Ne cannot be expected to be exceptionally 
well described by the RMF-I-RRPA theory. In few nu- 
cleon systems, the single particle structure around the 
Fermi surface is of the utmost importance in the calcu- 



lation of low-lying excitations. As such, the results pro- 
duced in a self-consistent mean field calculation are not so 
reliable. A better description would require a proper ac- 
count of excitations to the continuum above the coulomb 
barrier and probably for higher order correlations at the 
time dependent mean field level. The situation improves 
in heavier nuclei, were mean field theories were designed 
to yield good results at low computational costs. 
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TABLE I: ph structure for the 5.7 MeV Ml transition mode in 
^°Ne for the NL3. The second column refers to the normaliza- 
tion of the RPA amplitudes. The level quantum numbers in 
the third column are ±0.'^ , where ±f2 is the angular momen- 
tum projection over the symmetry axis and tt is the parity. 
In square brackets are the quantum numbers of the oscillator 
state which contributes most to the mean field single particle 
level. The effect of coulomb interaction can be seen as the 
small differences in the mixing percentages for protons and 
neutrons. A calculation with the electromagnetic interaction 
switched ofi' gives as a result a perfect isospin symmetry, with 
no differences observable within the accuracy of the computed 
results. 

Nevertheless, it is still interesting to delve further into 
the study of the properties of the main excitation peak, 
as the same analysis can be performed in heavier nuclei 
and many of the general features will still be present. 
The study of the structure of the excitation peaks can be 
carried out in detail attending to their ph structure. The 
contribution Cpu from a particular proton or neutron ph 
configuration to a RPA state is determined by 

c,h = {\x;,,\'-\y;:^?) (Ill) 

where X'^ and are the RRPA amplitudes associated 
with a particular excitation energy. Table [J outlines the 
single particle decomposition of the dominant Ml peak 
observed in Figure [T] All the strength is provided by 
a single particle transition within the sd-shell, from the 
last level in the Fermi sea to the first consecutive unoccu- 
pied level. The low collectivity indicates that, within the 
RMF+RRPA model, the spectrum of the Ml operator in 
^"Ne is of single particle character. Each of the two Dirac 
spinors in the particle-hole pair corresponds to a eigen 
state of the static RMF potential. They can be charac- 
terized by the Nilsson quantum numbers Vl'^[NnzK\ of 
their largest component in an expansion in anisotropic 
oscillator wave functions. Here O is total angular mo- 
mentum projection onto the symmetry axis, tt is the par- 
ity, N — 2nr + + A is the major oscillator quantum 
number, and A = — rus is the projection of the orbital 
angular momentum on to the symmetry axis. From these 
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quantum numbers one concludes the following approxi- 
mate selection rules: Afi — +1, AN — 0, Auz — —1 and 
A A = +1. The orbital character of the excitation peak 
is confirmed by the fact that Afl = AA, which implies 
that the change in the magnetic quantum number rus is 
zero. It is interesting to note that even if the approxi- 
mate selection rule points to an orbital character for the 
mode, the spin strength is nevertheless dominant. 
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FIG. 8; (Color online) Intrinsic transition density for neu- 
trons (left) and protons (right) of the Ml peak at 5.7 MeV 
in in ^"Ne. Full (red shade) and dashed lines (blue shade) 
indicate positive and negative values, respectively. The z- 
coordinate runs along the symmetry axis and r is the distance 
from the symmetry axis. The thin dotted line represents the 
rms radius of the ground state neutron or proton density, and 
qualitatively marks the position of the ground state nuclear 
surface. 

It is difficult to form a geometrical image of the nature 
of an excitation, having at hand only the information 
given in Table HI For that purpose, it is always useful to 
compare the the neutron and proton intrinsic transition 
densities in a plot. Figure [8] shows a color plot of the 
transition densities at an excitation energy of 5.7 MeV. 
Color is used to indicate the value of the function, with 
blue for negative values and red for positive ones. Re- 
gions with the same kind of line (solid or dashed), or 
color shade (red or blue), for both protons and neutrons 
are indicative of an in-phase vibration, while in regions 
where the opposite is true protons and neutrons vibrate 
out-of-phase. In this case the excitation is of clear isovec- 
tor nature, and we can observe the typical structure of 
a scissors mode; neutrons and protons are out of phase 
over the full space, with a concentration near the caps of 
the prolate nuclear shape. 

In such a simple case as the one found in •^"Ne the 
interpretation of the two dimensional color plot for the 
transition densities is very clear. They represent the in- 
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FIG. 9: (Color online) Radial part of the projected (to / = 
1, M = 1) transition densities of the Ml peak at 5.7 MeV. 
r is the distance from the symmetry axis. The prominent 
isovector nature is patent in the graphs. 



trinsic transition densities, referred to the intrinsic frame 
of reference, where only the total angular momentum pro- 
jection on the symmetry axis is well defined. In that re- 
gard, they are expected to contain admixtures from all 
possible angular momenta. However, the transition op- 
erator (Ml in this specific case) restricts the major con- 
tributions of the transition densities to the total response 
to its own total angular momentum, i.e. to in the case 
of Ml transitions to / = 1. It is therefore advisable to 
project out the weaker-contributing angular parts from 
the densities to obtain the actual transition density that 
would be observed in the laboratory frame of reference. 
For the Ml operator that means retaining , with the help 
of Eq. (|107p only the contributions coming from angular 
momentum 1 — \. In Figure [5] the radial part of such a 
projected transition density is plotted for the main peak 
in the ^'^Ne Ml response. 

Both transition densities are almost the mirror of each 
other, a very clear indication of the pure isovector nature 
of the mode at 5.7 MeV. We have already seen that in 
simple geometrical terms this mode can be interpreted 
as a rotation of neutrons against protons around an axis 
perpendicular of the symmetry axis. Furthermore, de- 
tails in Figure [S] show that two distinct regions can be 
distinguished. They are separated at around 2 fm from 
the origin, where the direction of rotation for protons 
and neutrons changes. The appearance of two regions 
(as depicted in Figure [§]) is already a strong hint that 
the simple picture of the proton density rotating against 
the neutron density as rigid rotors (as in the Two Rotor 
Model 64]) does not refiect reality in this nucleus. The 
traditional scissors picture considers the neutron and pro- 
ton densities as the blades of a scissor oscillating against 
each other. In addition, one has to take into account 
the angular momentum inherent in a ivT = 1+ excitation: 
it can be descibed as an oscillation of the scissor which 
rotates at the same time rotating slowly around its Ion- 
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gitudinal symmetry axis. However, the picture we derive 
from the results of our calculation is somewhat different. 



B. The Electric dipole (El) response 

We have chosen the electric dipole response in ^°Ne as 
a second example application of the RMF+RRPA for- 
malism with axial symmetry. In Fig. [TO] the El strength 
is plotted as calculated with the NL3 parameter set. The 
red curve corresponds to excitations along the symme- 
try axis with = 0^, while the blue curve are those 
perpendicular to the symmetry axis with K'^ = 1~. In 
principle, for prolate nuclei, as is the case for ^"Ne, the 
strength due to the K'^ = 0^ mode should lie at lower 
energies compared to the K'^ — \~ mode. As the nu- 
clear potential must be flatter (more extended) along the 
symmetry axis, it is energetically more favorable for the 
nucleons to oscillate in that direction than along an axis 
perpendicular to the symmetry axis, where the nuclear 
potential is narrower. It is possible, therefore, to relate 
the nuclear deformation with the energy separation of 
the two modes [80ll8ll|. 
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FIG. 10: (Color online) El strength in ^"Ne, as calculated 
with the NL3 parameter set. 

The splitting of the response due to the broken spher- 
ical symmetry, and its interpretation, can be observed 
in Figures [TT] and [TO] The former is the transition den- 
sity for the main IVGDR peak at 16.73 MeV observed 
in the K"" — 0~ response, while the latter corresponds 
to the peak at 21.31 MeV in the K"" = 1" mode. The 
prolate deformation is evident, as the intrinsic transition 
densities are elongated in the direction of the z-axis. The 
character of the K"^ = 1^ mode as a vibration along a 
perpendicular of the symmetry axis is easily recognizable 
in Figure 1121 By comparison, the transition density in 
Figure [TT] is then easily interpreted as a vibration along 
the symmetry axis. As it is expected for the IVGDR, the 
neutron-proton vibrations are out of phase over the same 
spatial regions. It is more evident still looking at their 
respective projections to the laboratory system of refer- 



ence, which are shown in the lower plots of Figures [TT] 
and [12] 
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FIG. 11: (Color online) ^''Ne IVGDR transition density for 
the = 0~ peak at 16.7 MeV, NL3 parameter set. 

Coming back to Figure [TO] the response in the energy 
region between 15 MeV and 25 MeV corresponds to the 
IVGDR. Its strength is heavily fragmented into several 
peaks in an energy interval of about 3-4 MeV for both ex- 
citation modes. The main contributions to the strength 
curve come from more than four different peaks. For ex- 
ample, the K"^ = 1~ IVGDR response is composed, be- 
sides the already mentioned peak at 21.3 MeV, by four 
additional major peaks, situated at 19.6, 20.2, 21.8, and 
22.4 MeV respectively. Their projected transition densi- 
ties, in Figure [TO] show that all of them can be classified 
as a vibration of the neutron density against the proton 
density. 

The decomposition in ph components of the main 
K'" = 0" and K"" ^ I- IVGDR peaks can be found 
in Table [TT] The characteristic AiV = 1 pattern of the 
IVGDR is present in both peaks. The high collectivity 
indicates a very coherent excitation pattern that fits into 
the properties of a giant resonance. This phenomenon 
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FIG. 12: (Color online) ^°Ne IVGDR transition density for 
the = 1~ peak at 21.3 MeV, NL3 parameter set. 
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FIG. 13: (Color online) Projected transition densities for ma- 
jor = 1~ peaks contributing to the electric dipole response 
in ^°Ne. 



can also be observed in the transition densities, where 
the coherent superposition of ph pairs is evident in the 
absence of wavefunction-like features, and is easily inter- 
preted in a macroscopic picture where the proton and 
neutron densities oscillate one against the other. The to- 
tal percentage of the classical TRK sum-rule exhausted 
between 10 MeV and 40 MeV for the calculated El re- 
sponse in ^°Ne is 111%. In a fully classical system the 
share of the strength exhausted by the — l~ mode 
should be double than that exhausted by the K'^ = 0~ 
mode; however, with 86% of the TRK sum-rule coming 
from the K"" = 1^ mode and 25% from the K'^ = 0~ 
mode, it is obvious that it is no longer true for quan- 
tum systems, even though we do not fully understand 
the mechanism behind this phenomenon. 

iC" = 0" peak at 16.73 MeV 
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peak at 21.31 MeV 
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TABLE II: Particle-hole structure of the two main IVGDR 
modes. N and P indicate a neutron or proton ph pair, respec- 
tively. The second column is the percentage of the contribu- 
tion of each particular ph excitation. The Nilsson quantum 
numbers, labeling the particle-hole, are shown in the next 
columns. The last column gives the energy of the excitation 
in MeV. 



VIII. CONCLUDING REMARKS 

In the present investigation we have formulated the 
relativistic random phase approximation (RRPA) on the 
basis of a relativistic mean-field (RMF) model having 
axial symmetry in a fully self-consistent way, i.e., the in- 
teractions used in both the RMF equations and in the 
matrix equation of the RRPA are derived from the same 
Lagrangian, i.e. the same energy functional. As it has 
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been shown, this self-consistency feature is of vital im- 
portance for the fulfillment of current conservation and 
the proper decoupling of spurious modes without further 
adjustments in the interaction. 

So far, pairing correlations have not been included. 
The inclusion of such correlations will allow the appli- 
cation of this method to a large number of investigations 
in medium and heavy nuclei, in particular, in a first step, 
for a systematic study of low-lying electric and magnetic 
dipole strength over large regions of the periodic table. 
Of course, the study of the nuclear response to other 
electric and magnetic multipoles is also open to scrutiny. 
Since the formulation of relativistic proton-neutron RPA, 
once the main building blocks presented in this docu- 
ment are present, is mostly trivial, its implementation 
opens the door to the wide area of nuclear spin-isospin 
excitations, in particular the Isobaric Analog Resonance 
(lAR) and the Gamov- Teller Resonance (GTR), but also 
for many types of weak processes such as the /3-decay and 
neutrino-reactions in deformed nuclei. 

In conclusion, the relativistic RPA formulated for ax- 
ially deformed systems represents a significant new the- 
oretical tool for a realistic description of excitation phe- 
nomena in large regions of the nuclear chart, which 
has been accessible so far only by relatively crude phe- 
nomenological models. Its development, and the sample 
application presented in this document, show that its fu- 
ture use in nuclear structure and astrophysics will pro- 
vide an valuable insight into very important, and still 
open, questions about the nature of nuclear interaction, 
collective response, deformation effects and cross sections 
relevant for astrophysical processes. 



Acknowledgments 

Helpful discussions with R. R. Hilton and D. Vretenar 
are gratefully acknowledged. P.R. thanks for the sup- 
port provided by the Ministerio de Educacion y Cien- 
cia, Spain. The paper has also been supported by 
the Bundesministerium fiir Bildung und Forschung, Ger- 
many under project 06 MT 246 and by the DFG clus- 
ter of excellence "Origin and Structure of the Universe" 
(www.universe-cluster.de) . 



APPENDIX A: TWO-BODY MATRIX 
ELEMENTS 

Starting form the general expression for the two-body 
matrix element in Eq. (|79p 



we obtain 



{kl'\Vt\k'l) = ' ^''^ 



(2^)^ 



{k\Q^{q)\k')^M{l\QM)\l' 



(Al) 

we first have to evalute the matrix elements ([50]) for the 
single particle opreators Q^{q). In cylindrical coordinates 

q = {qx,qy,qz) = {q cos x,q sin x,qz) (A2) 



Kk'iq,X,qz)^ I V:;,g„,r'^^fc,e^'?^^+^«"°^(^->^) 

(A3) 

It turns out to be useful to classify the various vertices 
by the spin quantum numbers S and S'^ = E of the ex- 
changed meson. For this reason we use the the 7-matrices 
in the Dirac basis defined by 

{7°, 7+ = ^(7' +*7'),7~ = 71^^' " ^^"^^ 



and obtain 

7^7m = 7"7" + 7 ' 7 +77'- 7"7" 



(A5) 



Including the scalar mesons (and neglecting for the mo- 
ment the isospin) we therefore have 5 vertices character- 
ized by the index fl: 



F'^ = (1,70,7+, 7-, 7^). 



(A6) 



fl runs over jx = s (for scaler mesons, 5 = E = 0), /i = 
(for the time-like part of the vector mesons, S — Yj — Q) 
and /i 3 (for the spatial parts of the vector mesons 

with S* = 1). The channels jl — ± describe the spin flip 
(E = ±1) and /i = 3 has E = 0. 

Using the Dirac spinors in cylindrical coordinates (j65p 
and exploiting the properties of the 7-matrices defined in 
Eq. (|A4p . we find that the (/^-dependence of the ampli- 
tudes i/ifc (r, Lp, z) -ij^ki (r, ip, z) can be separated 

V;fc(r,^,z)F'^^fe.(r,^,z) -z^Ff;^,(r,z)e'^^ (A7) 

where the integer 

A = Ofc - r^fc' - E = X - E (A8) 

is the orbital part of the angular momentum of the pair 
in z-direction. The real functions F^f,,{r^ z) are given by 



(All) 
(A12) 



This allows us to evaluate the iy9-integration in the inte- 
gral (|A3p analytically and to express it in terms of Bessel 
functions of the first kind 



Fkk' (^1 


^) 


= /^/^ 


+ fkfk'- 


-9t9t 


- 9k 9k' 


F'kk'ir. 


^) 


= f^f^' 


+ fkfk'- 


^9t9t' 


+ 9k 9k' 


Fkk'ir, 


^) 


= 9Uk^' 


- ft9l' 






Fkk' 


^) 


= fk9t' 


-9k ft' 






Fik'ir, 


^) 




-gift' 


^ fk9k' 


- 9k fk' 



We obtain 



tk'{q) = ^''^'e^''''J'L'{q^q^) 



(A14) 



(A15) 



with the integrals 



Hk'{q^q^) = / F^,,{r,z) JA(gr)e^«-. (A16) 
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which are either real or purely imaginary. Using the par- 
ity relation 

F4,(r,-z) = 7r(-)^+V^,,(r,z) (A17) 

whe recognize that the exponential factor e**^^ reduces 
to the cosine or sine depending of the quantum numbers 
TT and K of the mode and on the spin quantum numbers 
5, S of the vertex jl 



^^q.^^l cos(g.z) for 7r(-)^+s ^ = +1 . 
^ \isin(g,z) for tt{-)^+S~^ ^ -I ^^^^> 

Substitution of these expressions in the integral of 
Eq. (|Aip the ^-integration can be carried out analyti- 
cally and leads to the selection rule 



(A19) 



and to the following matrix elements for the exchange of 
scalar mesons 



kl'\VP''\k'l 



■ps A rs* 



(A20) 



of time-like part vector mesons 

{hl'\Vl^\m) = +J (A21) 
and of space-like of the vector mesons 

'kl'\VS'^\k'l) = J ^T+,A^Tu; (A22) 
'^'%T-„A^T+; (A23) 



(2^) 
(2^ 



2 kk'^^-' W 



3* 



.2 ■^kk''-^'^-^W 



(A24) 



where, for simplicity, we have neglected for each matrix 
element a factor Sq^. -n^, _n,/ and the arguments in the 
functions JF^^, {q, qz) and in the propagators 



A„i{q,qz) 



+ ql + 



(A25) 



in the two-dimensional momentum integrals. 

APPENDIX B: NON-LINEAR a PROPAGATOR 
IN MOMENTUM SPACE 

The equation to solve is 

[-A + ml + W{r)]5(j{r)^~gJp,{r) (Bl) 

with 

W{r) = 2g2a + Sg^a^ (B2) 



Because of axial symmetry and using cylindrical coordi- 
nates r = (r cosi^, r sin z) W{r) :— W{r,z) does not 
depend on the azimuth angle (p. We solve Eq. (jBip in 
momentum space. W{r) is local in r-space, but it is an 
operator in momentum space 

W{q,q')^ j (frW{r)e-''-^''-'i'^ (B3) 

The propagator in momentum space is the solution of 

(q2 + m^)A{q, q') -f j d\W{q, q")A{q", q') - 5{q - q') 

(B4) 

where the * is the convolution operator and W{q) is 
the Fourier transform of W{r). Expanding the 5- 
function in cylindrical coordinates in q-space using q — 
(q cos X, q sin x, qz) we find 



S{q-q') 



^-^^Siqz-q'z) E ^^"'"-"'^ 



(B5) 



n— — oo 



and taking the following ansatz for A 

oo 

A(g,g')= E ^n{q,q.,q',q'z)e"'^^-^'^ (B6) 



and inserting it in Eq. (jB4p leads to a set of integral 
equations for each A„ 

{q"^ + ql + ni'^)Xi{q,qz,q\qz) 

+ J (fq"Wniq,qz,q",qz)\iiq",qz,q',q'z) 

5{q-q') 



-Siqz - q'z), 



(B7) 



where we have used the obvious notation (Pq = qdqdqz- 
Each Wn can be calculated using the following series ex- 
pansion 



(B8) 



that leads to the following expression for the non-linear 
a field in momentum space 

Wn{q.qz,q',q'z)^ j ^W{r, z)e-^^'^'-'^> .Uqr)Uq'r) 



that together with Eq. IB 71 allow the numerical evaluation 
of the non-linear cr propagator in momentum space. 



APPENDIX C: EVALUATION OF THE Ml 
SINGLE PARTICLE MATRIX ELEMENTS 

The Ml operator is defined as: 
13 



fJ-N (gsS + gil) 



(CI) 
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with 



^ } protons 



9s 
91 



gn 



91 ^ ^ ' ' 3i = 

In the spherical coordinates defined as 



} neutrons (C2) 



x+ = -^{x + iy) , X- = ^{x- iy) , xq = z (C3) 



we find 



1 



2V2 



dr 



(C4) 
(C5) 



and the single particle matrix elements are 

(fc|£)o|fco) = eoff fc.,a,, j zd\ (D4) 

{fUk'+fkfk'+9t9i,+9f:9^d 
{k\D+\Q = eeff 5n,,n,,+i J rSr (D5) 

[fUi' + fifi' + 9t9l + 9t9l) 

where Ccff = eN/A for proton pairs and Ccfi = —eZ/A for 
neutron pairs. 



APPENDIX E: APPROXIMATE ANGULAR 
MOMENTUM PROJECTION 



and in the Afn single particle matrix elements the in- 
tegration over the azimuthal angle ip can be carried out 
analytically. This yields 



1 / 3 



- ^ {fidrfi, + ludrf^, + 9tdrgt, + g+drg^,) 

+ ^^^^ifU^'+9Ui) 



+ 



ifk fk' +9k 9k') 



{fUk-+9t9k') 



(C6) 



APPENDIX D: EVALUATION OF THE El 
SINGLE PARTICLE MATRIX ELEMENTS 

The effective isovector dipole operator, with spuri- 
ous translation of the center of mass already subtracted, 
reads in spherical coordinates 



p—l n—l 



(Dl) 



With the spherical coordinates of Eq. jCSj) the dipole 
operators are given in cylindrical coordinates as 



p—1 n—l 

N ^ 7 ^ 



(D2) 
(D3) 



The wave function \^im) in the laboratory frame is 
obtained by angular momentum projection from the in- 
trinsic wave function 1$) 



1^ 



IM) 



Y.9kPLk\'^) 



(El) 



K 



where the projector operator Plu^ is given by [s 



2/+ 1 

87r2 



dnvi}ji{n)R{n). (E2) 



r2 represents the Euler angles (ck,/3,7), Vj^jj^lfl) are the 
Wigner functions [s^l and R{n) 



_ p-iaJz p-il3Jy p~i^Jz 



IS 



the rotation operator. Taking into account the transfor- 
mation law for the multipole operators Qx^ under rota- 
tions 



(E3) 



The matrix element of this operator between two states 
with good angular momentum is given by 



with the reduced matrix element defined by 

OTT 



(E4) 



(E5) 



2^ (-) '9Kf9K. [ ^, ^ 



II, fj.' 



In the case of axial symmetry, the integral in the last 
equation is reduced to 



d(cos/3) d%,^^mKf\Qx^e~'P^y\K,) (E6) 
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To evaluate the overlap integrals in the last equation we 
use in the limit of the needle approximation (ssl . [s^ to 
first order in a Kamlah pH expansion 



-ipjy 



{Kf\Q^^\K,){Kf\e-^l'-^y\Kf) 

(E7) 

and using that the integral over (3 contributes only dX (3 — 
and TT we obtain the final expression for approximate 
angular momentum projection used in the calculation of 



RPA single-particle observables: 

{IfKf\\dx\\hK.)^{2U + l){2If + l) 



(E8) 
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